%% Load relevant data sets
load('522_epg_struct.mat');
load('boundaries.mat');

%% get info for recording3
% sampling rate for AUDIO
% 3RD acquisition
audioSR = new_epg_struct{3, 1}.Audio_sampling_rate;
epgSR = new_epg_struct{3,1}.EPG_sampling_rate;

audio = new_epg_struct{3, 1}.Audio;
%time = 0:(length(audio)-1) *  1/audioSR;
time = jshTime(audio,audioSR);
plot(time,audio)

% plot 1st two rows of sensors across time (1:14)
data = new_epg_struct{2,1}.EPG_data(:,1:14);

% number lit
sumEachRow = sum(data,2); % sum along columns (second dimension)
%equivalent to:  rowSums = data * ones(14,1);
size(data)
epgSR = new_epg_struct{1, 1}.EPG_sampling_rate;

mytime = jshTime(sumEachRow,epgSR);
plot(mytime,sumEachRow)

%% Write sound to WAV file and annotate in Praat
% wavplay.m doesn't work on Mac

% save all audio as WAV to annotate in Praat
wavwrite(new_epg_struct{2,1}.Audio,new_epg_struct{2,1}.Audio_sampling_rate,'audio2.wav')

for i=1:length(new_epg_struct)
   new_epg_struct{i}.Audio 
   wavname = strcat('epg_',num2str(i),'_',new_epg_struct{i}.File_name,'.wav');
   wavwrite(new_epg_struct{i}.Audio,new_epg_struct{i}.Audio_sampling_rate,wavname)
end

% THEN use Praat to annotate the boundaries...


% ... OR, the boundaries are actually given in the given data set,
% boundaries.mat, so in the interest of time, just use those.
% First column is a vector of vectors length 2 that represent the
% xmin and xmax of boundaries

%% Determine sample-start and sample stop numbers of regions of interest

% Convert boundaries to sample numbers
    
boundtime = cell2mat(boundaries);
boundsamps = round(boundtime*epgSR);
    % in boundsamps (4x4), rows correspond to struct number,
    % and columns in order: 
        % 1. palatal-startSample 
        % 2. palatal-stopSample
        % 3. [t]-start sample
        % 4. [t]-end sample


%% need function to convert sensor number to row,column indices 
% (and back again?)
sensors = [NaN,1:6,NaN,7:62];
sensmatrix = reshape(sensors,8,8)';


% % create 3-d matrix
% R = zeros([8,8,8]);  %repeat matrix
% 
% % 3rd dimension is so each row of sensors can be defined by an 8x8 matrix
% % where that row has ones while the others have zeros
% 
% %ids = nan(8,1);
% for i=1:8
%     R(i,:,i) = 1;
% %    idssensmatrix .* R(i,:,i);
% end


%% Calculate Contact Anteriority

% INVOLVES: Subset columns of samples of interest & sum rows

% Within each struct (trial?), compute the palate row sums of contacts that
% are activated for each sample

% OUTPUT is a matrix (8 x Nsamples) called epgRowSums for each struct;
% matrix row number corresponds to palate row number, where 1 is front, and
% 8 is back

for i=1:length(new_epg_struct)
   [Nsamples,Ncolumns] = size(new_epg_struct{i}.EPG_data);

   new_epg_struct{i}.epgRowSums = nan(8,Nsamples);
   new_epg_struct{i}.epgColSums = nan(Nsamples,8);
   for j=1:8
      % compute the RowSums
       new_epg_struct{i}.epgRowSums(j,:) ...
          = sum(new_epg_struct{i}.EPG_data(:,sensmatrix(j,isfinite(sensmatrix(j,:)))),2);
            % sum(X,dim) --> dimension to sum across (2 would sum across
            % columns, giving a sum for each row)

      % compute the ColSums
       new_epg_struct{i}.epgColSums(:,j) ...
          = sum(new_epg_struct{i}.EPG_data(:,sensmatrix(isfinite(sensmatrix(:,j)),j)),2);
      
   end
   
   
   % Compute Contact Anteriority for each sample all at once (acting on
   % whole matrix)
   new_epg_struct{i}.CA = contAnt(new_epg_struct{i}.epgRowSums);
   
   % Compute Centrality
   %    ...for left side
   new_epg_struct{i}.CC_L = contactCentral(new_epg_struct{i}.epgColSums(:,1:4));

   %    ...for right side
   new_epg_struct{i}.CC_R = contactCentral(new_epg_struct{i}.epgColSums(:,8:-1:5));

   
%    % alternative way to compute in a for loop
%    new_epg_struct{i}.CA = nan(1,Nsamples);
%    for k=1:Nsamples
%        new_epg_struct{i}.CA(1,k) = contAnt(new_epg_struct{i}.epgRowSums(:,k));
%    end
end

%% get Measures summary token-by-token (subsetting for interest times)

for i=1:length(new_epg_struct)
    sampSubset = boundsamps(i,:);
    new_epg_struct{i,1}.pal_CA = new_epg_struct{i,1}.CA(1,sampSubset(1):sampSubset(2));
    new_epg_struct{i,1}.t_CA = new_epg_struct{i,1}.CA(1,sampSubset(3):sampSubset(4));
    
    new_epg_struct{i,1}.pal_CC_L = new_epg_struct{i,1}.CC_L(sampSubset(1):sampSubset(2),1);
    new_epg_struct{i,1}.t_CC_L = new_epg_struct{i,1}.CC_L(sampSubset(3):sampSubset(4),1);
    
    new_epg_struct{i,1}.pal_CC_R = new_epg_struct{i,1}.CC_R(sampSubset(1):sampSubset(2),1);
    new_epg_struct{i,1}.t_CC_R = new_epg_struct{i,1}.CC_R(sampSubset(3):sampSubset(4),1);
    
end

%% Create EPG Graph

% take epg sample and turn into square matrix

mat1 = [0 0 new_epg_struct{1,1}.EPG_data(5,:)];
mat2 = reshape(mat1,8,8);
mat3 = ~mat2';
imshow(mat3)
